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We study the glassy super-rough phase of a class of solid-on-solid models with a disordered substrate 
in the limit of vanishing temperature by means of exact ground states, which we determine with a 
newly developed minimum cost flow algorithm. Results for the height-height correlation function 
are compared with analytical and numerical predictions. The domain wall energy of a boundary 
induced step grows logarithmically with system size, indicating the marginal stability of the ground 
state, and the fractal dimension of the step is estimated. The sensibility of the ground state with 
respect to infinitesimal variations of the quenched disorder is analyzed. 



PACS numbers: 64.60Fr, 64.70.Pf, 74.60. Ge 

There has been much interest recently in tha.pj;operties 
of crystal surfaces upon disordered substratesEfQ, in par- 
ticular the effect of the presence of pinning potentials on 
surface roughening. Further motivation comes from the 
close relationship of the corresponding disordered solid- 
on-solid model to the two-dimensional vortex free XY 
model with randorp-fields and to the randomly pinned 
planar flux arrayETtj. The latter is particularly relevant 
with respect to the technologically important aspect of 
flux pinning in high-T c superconductorscJ. 

It has been demonstrated that these systems have a 
phase transition at a critical temperature from a ther- 
mally rough phase into a super-rough phase at low tem- 
peratures. Whereas the existence of this transition is 
established by now, the qualitative and quantitative fea- 
tures of the glassy low temperature phase are still de- 
bated. Predietiojis of earlier renormalization group (RG) 
calculationsc3i§Elta turned out to be inconaaatible with re- 
sults of extensive numerical simulations03trlLil. The sub- 
sequent discovery of the relevance CiLreplica sjaniactry 
breaking (RSB) effects in variationaO and RGE3E3 cal- 
culations lead to a variety of new results, from which 
those obtained by the variational treatment are agaiO|i|i 
disagreement with the most recent numerical studyU'lla. 

For us the situation seems to be the following: close 
to the transition the discrimination between various pre- 
dictions is numerically hard because of its smallness on 
intermediate length scales. Far from the transition (i.e. 
deep in the glassy phase at low temperatures) the pic- 
ture should become much clearer, meaning that the dis- 
order dominated effects become stronger. The clear- 
est evidence for the latter can be expected for strictly 
zero temperature, when the roughness is exclusively pro- 
duced by the substrate alone. However, to reach this 
limit Monte Carlo simulations of glassy systems like those 
we are interested in suffer from notorious equilibration 
problemalj. At zero temperature the properties of the 
system under discussion are completely determined by 
its minimal energy configuration or ground state. Since, 



as we show in this paper, it is possible to calculate this 
state exactly, a detailed and reliable picture of the zero 
temperature limit of the glassy phase can be obtained. 

Thus the aim of the present paper is twofold: first 
we calculate numerically the zero temperature limit of 
the height-height correlation function for finite systems 
and compare the result with various analytical predic- 
tions and finite temperature simulations. Second, mo- 
tivated by the observation of manifestly glassy features 
(like RSB, slow dynamics, metastable states) in the low 
temperature phase, we ask how far concepts developed 
for finite-dimensional spin glasses can be applied here 
and explore the nature of this ground state in much 
greater detail. We investigate its stability with respect to 
step-excitations, analyze this fractal boundary induced 
domain wall itself and study the chaotic nature of the 
ground state by an application of infinitesimal variations 
of the substrate heights. These important issues have 
never been discussed in the present context. 

The solid-on-solid (SOS)-model we consider here is de- 
fined by the following Hamiltonian: 

H = ^2f(hi-hj) (1) 

m 

where (ij) are nearest neighbor pairs on a d-dimensional 
lattice (d = 1,2) and f(x) is an arbitrary convex 
if"( x ) > 0) and symmetric (f(x) = f(—x)) function , 
for instance f(x) = x 2 . Each height variable hi = di +n, 
is the sum of an integer particle number, which can also 
be negative, and a substrate offset di £ [0, 1[. For a flat 
substrate, (L.— Q for all sites i, we have the well known 
SOS-modelE3. The disordered substrate is modeled by 
random offsets di £ [0, lp. In the present paper we are 
only interested in uncorrelated disorder, meaning that 
all offsets are distributed independently. The method 
we use is, however, applicable to any disorder distribu- 
tion, in particular to the case of correlated disorder (i.e. 
{didi+ r ] av — g(r), with g(r) an arbitrary function), toco. 
In what follows [. . .] av denotes the disorder average. 
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The random offsets induce a frustration into the system 
in the same sense as quenched randomness in the interaa- 
tion strengths does in the context of spin glass modelslia. 
The minimum energy surface (being the set of particle 
numbers {n{\ that minimize the energy function H), is 
no longer flat (hi = const.) but a highly non-trivial ob- 
ject and the calculation of this ground state is a complex 
combinatorial optimization problem. By introducing the 
height difference variables for links on the 

dual lattice we can reformulate our task as a minimum 
cost flow problem 



Minimize c 

(ij) 



*? ( x ij ) 



(2) 



with intcg and convex, flow dependent (i.e. in- 

dependent) cost functions Cij(x) = f(x — dij), dij = 
di — dj. For this problem we developed an efficient 
pseudo-nolynomial algorithm which is described in much 
detail irEd. It is guaranteed to find the exact ground 
state, typically in 15 minutes on a Sparc 20 workstation 
for a system with N — 256 x 256 sites. We used fixed 
boundary conditions (b.c.) for the model ([l] becuase for 
technical reasonsEd the algorithm works only on a planar 
graph. We would like to point out that this is the first 
extensive application of a minimum cost flow algorithm 
to a problem in theoretical physics to our knowledgeEj. 

Much effort, numerically as well analytically, has been 
devoted to the calculation of the height-height correlation 



function C(r,r') = [(^nnf 
close to the transitionErHr^ 1 



n the case f(x) 



The zero tempera- 
ture limit has not been investigated so far and the results 
we present here are the first reported in the literatureo. 
At high temperatures (T ;» T g = 2/n) one expects 
C(r, r') — T/27rlog |r — r'| in an infinite lattice. Because 
of the natural bending of this logarithmic curve by any. 
boundary conditions this expression should be replacedtll 
for finite lattices by C(r,r') = T/2TrPi,(r,r'), where for 
fixed b.c. the propagator is given by 



P L (r,r') 



2 
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E 



(3) 



[sm(qoxn) sm(qoym) — sin(qox'n) sin(qoy'm-)] 2 
2 — cos qon — cos q^m 

with r = (x, y). r' = (x 1 , y') and qo = ir/(L + 1). As men- 
tioned above at low temperatures T < T g different sce- 
narios are currently under discussion in the literature: a 
Gaussian scaling, with C(r) linear in log(r), and a log 2 (r) 
behavior. In the first case a plot C(r) versus Pl{t) should 
yield a straight line, in the second case one should be-able 
to fit C(r) to a quadratic polynomial in Pl{t) (seetZl). 

We define the site averaged correlation function: 
C(r) = 2/ L 2 J2t=iJ2y=i[( h {*,y) ~ h (x+r, y ))\-v and also 
the corresponding site averaged lattice propagator Ph{r), 
which behaves for L — > oo and r <C L like Pl{t) ~ 



log(r)/27r. In fig. 1 we show C(r) versus Pl{t) for a sys- 
tem size L = 128. Obviously one does not get a straight 
line, which one would obtain if C(r) ~ oq + a\Pi,(r). 
Thus the Gaussian scaling C(r) oc log(r) for r — > oo can 
definitely be excluded on these length scales. 
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FIG. 1. The site averaged correlation function C(r) ver- 
sus the lattice propagator Pl{t) for L = 128 and aver- 
aged over 2000 samples. The broken line is a least square 
fit to C(r) = 0.008 + 0.2LP £ (r) + 0.57P L (r) 2 . The inset 
shows C(r)/Ph{r) versus Pl(t), the straight line indicates 
the amount of curvature of the data. 

Moreover, although a fit like C(r) ~ a + aiPi(r) + 
ci2PL{r) 2 works fine on first sight, there is still a signif- 
icant bending in the plot of C(r)/ Pl{t) versus Pl{t) 
shown in the inset of fig. 1, indicating a possible pres- 
ence of even higher order terms. The fit parameters 
ai = 0.21 and a-i = 0.57 are compatible with a lin- 
ear zero temperature ex-fcsapolation of the finite T re- 
sults of Marinari et alO. The. -existence of a log 2 - 
term is compatible with earlicraEl'aEl replica symmetric 
RG esdeulations and with the most recent one involving 
RSBll3ll3. Note that the RG prediction for the coefficient 
isEI 02/ 2,11 = 2/tt 2 ■ t 2 + 0(t 3 ), where r is the distance 
from the critical point. Since r = 1 in our case (T = 0) 
it is hard to compare the numerical value of the coeffi- 
cient. Neglecting the higher order corrections one gets 
0,2 /27r = 2/-7T 2 k 0.20 which is twice as large as our esti- 
mate. 

We also studied the correlation function for different 
convex cost functions f(x) = \x\ n with varying n. We 
observed that with increasing power n for smaller dis- 
tances large height differences \xtj\ > 1 are suppressed 
due to their larger costs. At larger distances, however, 
the roughness increases systematically with n, expressed 
in monotonically increasing fit values for the coefficient of 
the log 2 -term. This is due to the decreased costs for small 
height differences (xij = ±1). Moreover, we conclude 
that the coefficient for the log 2 -term is a non-universal 
number at zero temperature because of its significant de- 
pendence on the actual shape of the cost function f(x). 

An intriguing question concerns the stability of the 
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ground state with respect to thermal fluctuations. To 
attack this problem we take over the concept of domain 
wall renormalization that-is well known in the context of 
random spin systemsEjcj: we ask how much energy a 
step (or domain wall) of height one in a system of lin- 
ear size L would cost. If this energy is an increasing 
function of the size L, it indicates the stability of the 
ground state even at finite temperatures (disregarding 
other, more complicated excitations). The step is in- 
duced by appropriate boundary conditions: we fix the 
lower boundary to zero (fy^o) — 0) an d the upper bound- 
ary to one (/i^l+i) = !)■ At the left and right boundary 
we enforce fyo.t,) = h(L+i,y) = Q(y — L/2), with 6(x) = 
for x < and 0(x) = 1 for x > 0. This procedure in- 
duces a straight step {ht x ,y\ = for y < L/2, hi x , v ) = 1 
for y > L/2) in the pure case (di = 0), which costs an 
energy AE — L (implying that the flat surface is sta- 
ble at finite temperatures in this case). In addition to 
quadratic geometries we also varied the height (i.e. the 
distance between the upper and lower boundaries) and 
considered rectangular geometries with L x H sites. 

We calculated the ground state first for hi — on the 
whole boundary and then for the step inducing bound- 
aries described above (we choose fix) = x 2 from now 
on). In fig. 2 we show the averaged step energy as a 
function of system size. The data are nicely fitted by a 
logarithmic L-dependence 
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61nL with b = 0.52 ± 0.02 



(4) 



At first sight the tempting conclusion of this observation 
would be that there exists a low temperature phase with 
long range order, with an order parameter given by the 
overlap of the system with its ground state (note that the 
boundaries are fixed so that this definition is unambigu- 
ous). However, due to the strong chaotic rearrangement 
of the state upon temperature changes, which we discuss 
below, most probably this quasi long range order gets de- 
stroyed at finite temperatures. We think that this result 
merits a further investigation e.g. via simulations. 

An inspection of the geometrical shape of the step that 
is induced by the boundary condition described above re- 
veals that it is fractal, crossing the whole sample. There- 
fore we varied H to ensure that the limitation in the 
y-direction does not influence the quantitative results we 
obtain. We calculated the averaged length of the step as 
a function of system size, which is depicted in the inset 
of fig. 2, and obtain a good fit to given by 

£stc P oc L dstep with 4tcp = 1-35 ± 0.02 . (5) 

where d s t op is the fractal dimension!! of the step. This 
is the first estimate for such an exponent in the present 
context, which we expect to be universal, as we checked 
for various forms for f{x). 

In order to shed further light on the glassy features 
of the low temperature phase of our system we study 
what became known under the notion of chaos in spin 
glassesESO. For the latter one observes that infinitesimal 
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FIG. 2. The averaged step energy [|Ai5|] av as a func- 
tion of system size L, averaged over 10 4 samples. It is 
fi(L) = 0.15+0.56 ln(L) and f 2 (L) = 0.19+0.52 ln(L), see eq. 
(^). The inset shows the step length C s tc P as a function of sys- 
tem size. Here it is gi (L) = 0.77L 1 ' 37 and g 2 (L) = 1.04L 1 ' 33 , 
see (||). The boxes are data for an L x 2L-geometry, the 
crosses for L x L. Note that for larger H the step energy is 
slightly smaller and the step length slighly larger, since the 
step has more space to optimize its configuration. 



variations of temperature or interaction strengths decor- 
relate the state of the-system over a distance that is called 
the overlap lengthESEJ. Here we study the same scenario 
by infinitesimal (random) perturbations of the offsets of 
the sample and comparing the resulting changes in the 
ground state configuration. To be concrete we define new 
offsets by d\ = di +£j, with £j 6 [—8/2, +6/2} and S <C 1 
and given the unperturbed and perturbed ground state 
heights hi and h! { , respectively, we are interested in the 
accumulated height difference 
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To get an idea about the scaling behavior of this quan- 
tity let us consider the one-dimensional case first, for 
which the ground state (with free right boundary con- 
dition) can be constructed iteratively. The iteration 
/ij+i = H(hi, di+\ — di) leads asymptotically to a ran- 
dom walk: [{hi-\- r — /ii) 2 ] a v ~ r with a slightly modified 
short distance behavior. The same holds for the random 
variable hi — h! i (with the site index i as "time" ) , but its 
amplitude is reduced by factor S. Thus in one dimen- 
sion we expect Xi^(^) ~ 8 2 L, which we verified also 
numerically. Note that S is the inverse length scale be- 
yond which hi — h! i becomes typically different from zero 
if h = h' = 0. 

Motivated by this observation we hypothesize the fol- 
lowing scaling form for xl- 



Xl(5)/S^ ~ g(L) 



(7) 



with rj = 2 in one dimension and g(L) an arbitrary func- 
tion. In fig. 3 we show a scaling plot of the data we ob- 
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tained from our our ground state calculation. The data 
collapse is acceptable for ?y = 0.91 ± 0.05. One should in- 
terpret 1 /8 V as the characteristic length scale over which 
the two ground states decorrelated (in analpgs£-.to the 
overlap length in the context of spin glassescS'Ej) . Wc 
note that this is the first quantitative description of chaos 
is this disordered SOS model. 
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FIG. 3. Scaling plot of Xl as defined in (|6|) for various 
values for the perturbation amplitude 3, averaged over 10 4 
samples. The best data collapse is obtained for r\ = 0.91. The 
full line represents the scaling function involving a dominant 
log 2 L-term. 

To summarize, our numerical investigation of the zero 
temperature limit of the glassy phase in SOS models on 
disordered substrates revealed the following picture: 1.) 
There is a dominant log 2 (r)-contribution in the height- 
height correlation function as predicted by RG calcula- 
tions. However, we find also indications for the existence 
of higher nonlinearities in log(r). 2.) Different forms 
for the microscopic interactions yield qualitatively simi- 
lar but quantitatively different results for e.g. the coeffi- 
cient of the log 2 term (at T = 0). 3.) the step excitation 
energy increases logarithmically with the system size in- 
dicating a marginal stability of the ground state against 
thermal fluctuations. 4.) The boundary induced step it- 
self is fractal with a fractal dimension of df ~ 1.35. 5.) 
The concept of the chaos known from spin glasses also 
applies here, the chaos exponent for the overlap length is 
close to one. 

As a future perspective we would like to remark that 
with our minimum cost flow algorithm it is for the 
first time possible to attack efficiently a large variety of 
combinatorial optimization problems occuring for disor- 
dered and frustrated systems in statistical and solid state 
physics, in particulat-those involving "rea.1" flows like en- 
sembles of flux linestfl in lattice modelsE3. 
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